
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <time.h>
#include <omp.h>
//////////////////////////////
#include <cuda.h>
#include <cuda_runtime.h>
//////////////////////////////


#define REPEAT 256

__inline__ __device__ float warpReduceSum(float val) {
  for (int offset = 16; offset > 0; offset >>= 1) 
    val += __shfl_down_sync(0xffffffff, val, offset);
  return val;
}

__inline__ __device__ float blockReduceSum(float val, int tid) {
  static __shared__ float temp[32]; // Shared mem for 32 partial sums (32 = max warps)
  int lane = tid % 32;              // thread id within the warp
  int  wid = tid / 32;              // warp id

  val = warpReduceSum(val);         // Each warp performs partial reduction
  if (lane==0) temp[wid]=val;       // thread 0 of each warp writes warp-level sum
  __syncthreads();                  // Wait for all partial reductions to complete

  if (wid==0) {                     // Done only by warp 0
    val = (tid < 1024 / 32) ? temp[lane] : 0.0f;
    val = warpReduceSum(val);       // Call again warp-reduction to sum over the blocks
  }

  return val;
}

__global__ void grid_sum_naive(const float *A,float *sum){
  int id    = blockIdx.x * blockDim.x + threadIdx.x;
  atomicAdd(sum,A[id]);
}

__global__ void grid_sum(const float *A,float *sum){
  int id    = blockIdx.x * blockDim.x + threadIdx.x;
  float tmp = A[id];
  __syncwarp();
  tmp = blockReduceSum(tmp,threadIdx.x);
  if (threadIdx.x==0) atomicAdd(sum,tmp);
}


int main(void){
  long unsigned int N = 1024L*1024*512;
  float *a = (float *)malloc(N*sizeof(float));
  struct timespec ts0,ts1;
  double time;
  float sum;
  float *d_sum;
  
  fprintf(stderr,"\nFilling array\n");
  for (long unsigned int i=0;i<N;i++){
    a[i] = (float)(rand()-RAND_MAX/2)/(float)(RAND_MAX);
  }

  
  // Allocate arrays on GPU
  float *d_a;
  cudaMalloc(&d_a,N*sizeof(float));
  cudaMalloc(&d_sum,sizeof(float));
  
  // Transfer initial data to GPU
  cudaMemcpy(d_a,a,N*sizeof(float),cudaMemcpyHostToDevice);

  
  
  clock_gettime(CLOCK_REALTIME,&ts0);
  fprintf(stderr,"\treduction on HOST");
  for (int r=0;r<REPEAT;r++){
    sum = 0.0f;
#pragma omp parallel for num_threads(6) reduction(+:sum)
    for (int i=0;i<N;i++) sum += a[i];
  }
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-ts0.tv_sec + (ts1.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr,"       : %6.2f seconds, sum=%12.6f\n",time,sum);
  

#define B 256

  
  
  clock_gettime(CLOCK_REALTIME,&ts0);
  fprintf(stderr,"\tnaive reduction on GPU ");
  for (int r=0;r<REPEAT;r++){ 
    cudaMemset(d_sum,0,sizeof(float));
    grid_sum_naive<<<N/B,B>>>(d_a,d_sum);
  }
  cudaMemcpy(&sum,d_sum,sizeof(float),cudaMemcpyDeviceToHost);
  cudaDeviceSynchronize();
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-ts0.tv_sec + (ts1.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr," : %6.2f seconds, sum=%12.6f\n",time,sum);

  
  
  clock_gettime(CLOCK_REALTIME,&ts0);
  fprintf(stderr,"\treduction on GPU ");
  for (int r=0;r<REPEAT;r++){ 
    cudaMemset(d_sum,0,sizeof(float));
    grid_sum<<<N/B,B>>>(d_a,d_sum);
  }
  cudaMemcpy(&sum,d_sum,sizeof(float),cudaMemcpyDeviceToHost);
  cudaDeviceSynchronize();
  clock_gettime(CLOCK_REALTIME,&ts1);
  time = (ts1.tv_sec-ts0.tv_sec + (ts1.tv_nsec-ts0.tv_nsec)/1e9);
  fprintf(stderr,"       : %6.2f seconds, sum=%12.6f\n",time,sum);
 
  exit(0);
}
